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Abstract. We have performed 3-D numerical magneto- 
hydrodynamic (MHD) jet experiments to study the insta- 
bilities associated with strongly toroidal magnetic fields. 
In contemporary jet theories, a highly wound up mag- 
netic field is a crucial ingredient for collimation of the 
flow. If such magnetic configurations are as unstable as 
found in the laboratory and by analytical estimates, our 
understanding of MHD jet driving and collimation has to 
be revised. 

A perfectly conducting Keplerian disc with fixed den- 
sity, rotational velocity and pressure is used as a lower 
boundary for the jet. Initially, the corona above the disc 
is at rest, permeated by a uniform magnetic field, and is 
in hyd rostatic equilibrium in a softened gravitational field 



from a point mass. The mass ejection from the disc is sub- 
sequently allowed to evolve according to deviations from 
the initial pressure equilibrium between disc and corona. 

The energy equation is solved, with the inclusion of 
self-consistently computed heating by viscous and mag- 
netic dissipation. We find that magnetic dissipation may 
have profound effects on the jet flow as: 1) it turns on in 
highly wound up magnetic field regions and helps to pre- 
vent critical kink situations; 2) it influences jet dynam- 
ics by re-organizing the magnetic field structure and in- 
creasing thermal pressure in the jet; and 3) it influences 
mass loading by increasing temperature and pressure at 
the base of the jet. 

The resulting jets evolve into time-dependent, non- 
axisymmetric configurations, but we find only minor dis- 
ruption of the jets by e.g. the kink instability. 
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1. Introduction 

Though astrophysical jets have been observed in a rather 
wide range of accreting systems, it is generally assumed 
that the mechanism for acceleration and collimation is 



generic (Livio, 1997; Spruit, 2000, for recent reviews) 



Contemporary jet theories rely on magnetic forces as the 
jet producing mechanism, either with the flow emanating 
from the disc surface and associated with an op en mag- 
netic field structure ( Blandford and Payne, 1982 ) or with 
the flow emanating from the disc-star interface and asso- 
ciated with a closed field structure connecting disc and 
central star ( Pringlc, 198S ). In the former so-called disc- 
wind scenario, which is currently receiving most attention, 
the driving force, acting to overcome the gravitational pull 
of the central object, is typically regarded as a component 
of the centrifugal force along the magnetic field. However, 
adopting an incrtial frame of reference the driving mech- 



anism may equally be interpreted as magnetic (Spruit 



1996). The inertia of the outflowing gas eventually be- 



comes dynamically important and the gas stops corotat- 
ing with the underlying disc. At this point the magnetic 
field, which cannot easily slip through the highly ionized 
gas, is wound up in a cork screw manner between the 
disc and the vertically accelerated and less quickly rotat- 
ing outflow. The tension ("hoop stress") of the wound-up 
magnetic field configuration is generally believed to be the 
collimating force, but collimation by a poloidal field has 



been proposed as well (3pruit ct al., 1997) motivated by 
instability arguments against the wound-up field scenario. 

The assumed large scale magnetic field in the disc- wind 
scenario may either be provided by dynamo processes in 



the disc (Brandenburg, 200C; Bardou et al., 2001) or cap- 



tured from the environment and advected inwards by the 



accreting matter (Cao and Spruit, 1994). Numerical work 



has been carried out to investigate situations in which a 
signi ficant fraction of the fi e ld lines loop back into the 



disc ( Romanova et al., 1998 ; Turner et al., 199E ). These 



experiments have focused on the ejection and acceleration 
mechanisms as well as the evolution of the magnetosphere 
close to the disc, and have not followed the flow further 
out for extended times. 

The numerical work done in relation to the "classi- 
cal" large scale open field scenario may be divided into 
experiments attempting to include (par ts of) the accre- 
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tion disc in the computational domain (Uchida and Shi 



bata, 1985; Bell and Lucek, 1995) and more recent work 
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(Ustyugova et al., 1995; Ouyed et al., 1997; Meier et al. 
1997j) ^hich has used the disc only as a fixed lower bound- 
ary to avoid problems with radial collapse of the disc and 
thereby relatively short time evolution of the experiment. 
Both types of experiments have been axisymmetric and as 
such have not questioned the potentially (kink) unstable 
wound up magnetic field configuration on which the colli- 



mation process relies. However, Spruit et al. (2001) have 



recently proposed fast reconnection processes in more dis- 
ordered non-axisymmetric large scale magnetic fields in 
the jet like gamma ray burst (GRB) scenario. The mag- 
netic dissipation is proposed for the GRB fireballs e.g. to 
produce the observed radiation with better efficiency. We 
find that such magnetic processes may have severe impact 
on jet dynamics and stability in particular. 

The main purpose of the work presented here is to es- 
tablish whether or not the jets in the disc-wind scenario 
are prone to catastrophic MHD instabilities. This calls for 
an implementation of a setup in three dimensions which 
will be described in Sect. ||. For simplicity, we assume a 
large-scale open magnetic field structure and use the disc 
as a fixed boundary. The details of the initial conditions 



and boundary conditions are found in Sect. |2.2| and 2.5 re- 
spectively. Results are presented in Sect. [|, with Sect. |3~l| 
concentrating on the jet dynamics and Sect. 3.2 on the ob- 
served 3-D jet stability in the experiments. In Sect. ^ we 
discuss the results with special emphasis on thermal prop- 
erties and magnetic field structure. Finally, the paper is 
summarized and conclusions are presented in Sect. |^. 

2. Model 

The jet flow is described numerically by solving the MHD 
equations in the following form; 



dp „ 



(1) 



d(pu) 



Of 

fiJ = V x B, 



= —V • (puu — t') + VP + J x B - pV$, (2) 

(3) 



E = i]J — u x B, 



dB 

~dt 

de 
di ' 



= -V x E, 



-V-eu-PV-u + Q, 



(4) 
(5) 

(6) 



where p is the mass density, u the velocity vector, t' the 
viscous stress tensor, P the thermal gas pressure, $ the 
gravitational potential, B the magnetic flux density vec- 
tor, J the electric current density vector, E the intensity 
of the electric field, rj the electric resistivity, e the inter- 
nal energy per unit volume and Q the sum of viscous and 
Joule dissipation. 



2.1. Dimensions and numerics 

Both for numerical reasons and to be able to adapt the cal- 
culations to stellar as well as galactic scales, the physical 
quantities are given in units of characteristic quantities of 
the system. The form of the equations are not changed by 
converting to this system of units. The appropriate units 
may be constructed from the mass of the central object, 
M, the magnetic flux density, B, and the length scale, 
A. The characteristic length scale A is assumed to corre- 
spond to the inner disc radius. If not otherwise stated, the 
dimensionless quantities listed in Table |l| are used in the 
following sections. 

The code uses a sixth order accurate method for par- 
tial derivatives and a fifth order accurate interpolation 
method with the variables represented in non-uniform 



(Eulerian) staggered meshes (Nordlund and Galsgaard 
|l997a; Rdgnvaldsson, 199S ) . In the discrete representa- 



tion of quantities, the scalar variables (e and p) are zone 
centered and components of the vector variables (p = pu 
and B) are face centered in a unit mesh volume. The solu- 
tion for the eight variables, p, e, p x , p y , p z , B x , B y , B z is 
advanced in time by a third order predictor-corrector pro- 



cedure (Hyman, 1979), modified to accommodate variable 
time steps. The code has previously been verified by a 
number of standard test problems and hencef orth been 
used for exp eriments involving 3-D turbulence ( Nordlund 
|ct al., 19*94 ), investigations of problems rel ated to coronal 
heating ( Nordlund and Galsgaard, 1997b), buoyan t rise 
of magnetic flux tu bes (Dorch and Nordlund, 1998 ), dy- 
namo experiments (Dorch et al., 1999), stellar convection 
( Nordlund and Stein, 2000 ) and magnetized cooling flows 
( Rdgnvaldsson et al., 2002j ). 

Cartesian coordinates (x, y, z,) are used to describe the 
disc-corona system in the computational box. To calculate 
the derivative or the interpolated value at a given point, 
the six nearest grid points are involved. As exactly the 
same operators are used in the whole grid, three layers 
of ghost zones must be added to prevent periodic wrap- 
ping of the computational domain. For (obsolete) rea- 
sons of computational efficiency, we use x (first index) 
as the vertical direction and place ghost zones in the in- 
dex range i € [1 : 3] A [m x — 3 : m x ]. The y- and z- 
directions are taken to be periodic. The (y, z) cross sec- 
tion of the computational domain is quadratic, both in 
terms of number of grid points [m y — m z ) and physical 
size (L y — L z ), and centered so that y £ [—L y /2,L y /2] 
and z G [— L z /2, L z /2\. Only odd numbers of m y and m z 
have been chosen, making (x, y, z) = (0, 0, 0) at zone cen- 
ter (i,j,k) = (4, (m y + l)/2, (m z + l)/2). 



2.2. Initial conditions 

To investigate the magnetic driving and collimation of out- 
flows it is desirable to choose an initial state in which 
all other dynamical effects are small. Therefore, to obtain 
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Quantity 


Unit 


YSO 


AGN 


Velocity 
Time 

Mass density 


\J GM/X 
v/A3/(GM) 
B 2 \/(fj, GM) 


138 km/s 
0.58 days 
4.2 x 10" 14 g/cm 3 


0.2 c 
0.51 days 
1.8 x 10" 17 g/cm 3 



Table 1. The magnetic permeability is denoted by fio- Representative values are listed for a young stellar object (YSO) 
with M = IMq, B = 10 G, and A = 10 -Rq. Likewise, an entry is included for a typical active galactic nucleus (AGN) of 
M = 10 8 Mq, B = 100 G and A = 10 Rs, where the Schwarzschild radius is given by Rs = 2GM/c 2 . 



such a state, the corona is assumed at rest and in hy- 1995 ). Due to the reduced rotation velocity obtained when 



drostatic equilibrium in the gravitational field of a point 
mass ( |Ouyed and Pudritz, 1997 ). Self-gravity of the gas is 
neglected. A current free configuration is chosen so that 
the Lorentz force vanishes. In this way, the momentum 
equation for the corona reduces to 



VP = -pV$. 



(7) 



By the use of a polytropic equation of state, P = Kp 1 , the 
reduced momentum equation may be integrated, giving 
the density distribution for the corona, 



P = 



7^1. 



(8) 



A positive value of the integration constant, <&o > 0, may 
be introduced to keep finite densities at large distances. 
The value of the polytropic exponent is chosen to be 7 = 
5/3, corresponding to an adiabatic mono-atomic ideal gas. 

For numerical reasons, the gravitational potential is 
smoothed by introducing a so-called softening parameter, 



1 



V 



(9) 



Due to the location of the gravitational potential in the 
staggered meshes, a non-zero softening parameter is neces- 
sary to make the potential well-defined at x = y = z = 0. 
Furthermore, a softened potential makes the potential gra- 
dient and the initial density gradient less steep close to the 
central object and therefore easier to handle by the nu- 
merical code. The softening parameter was chosen to be 
l s = iy/l/2, for the experiments actually reported herein. 

The disc surface (lower boundary) is assumed to have 
a density distribution described by Eq. ||, and the corona 
is in pressure balance with the underlying disc. In the case 
of negligible thermal pressure, the softened potential gives 
rise to an equilibrium Keplerian like velocity structure in 
the disc given by 



(y 2 + z 2 + Z s 2 )3/4 ' 

y 

(y 2 + z 2 + Z 2 ) 3 / 4 



(10) 



Numerical experiments relying on a smoothed or softened 
potential produce generally weaker jets ( Bell and Lucck 



softening is introduced a longer time scale for the build 
up of the azimuthal field (and the jet flow) is anticipated. 
Devising identical experiments and only varying soften- 
ing, we find no significant effect in jet features but only 
a scaling of dynamical quantities. The main effect of soft- 
ening is an increase of the rotational period which scales 
the characteristic time and the specific kinetic energy (cf. 
Fig-0). 

In the chosen numerical setup, a density jump at the 
disc-corona interface must be applied with some caution, 
as the interpolation of the zone centered values could 
artificially bring mass into the corona. A density jump 
may however be implemented to have effect on the en- 
ergy flux as is demonstrated in Sect. |2~J| below. The pres- 
sure balance between disc and corona together with the 
assumed disc density structure expressed by Eq. ^ im- 
plies that a density jump may equivalently be regarded 
as a jump in internal energy per unit mass, T ~ P/p- 
Hence, the jump £ = (T COIona /T^ isc ) Q at the disc-corona 
interface expresses a jump in (isothermal) sound speed as 

Wel1 ; C s,disc = c s,corona/C- For C 2 ,di S c < K + U l) '> the ra " 

dial pressure gradient in the disc, dp/dr ~ (pc 2 )di sc /-R, 
becomes negligible compared to the centrifugal force and 
the pressure gradient does not influence the radial balance 
of the disc. Anyway, since the disc is not simulated, the 
disc structure will not be of a numerical concern. 

A constant vertical magnetic field, B x = Bq, is chosen 
to penetrate both disc and corona initially. In a more re- 
alistic case, one may expect the field to be inclined with 
respect to the disc surface and the field strength is likely 
to increase towards the center due to advection by the 



accreting matter (Cao and Spruit, 1994). Because of the 



periodic boundaries, such a magnetic field configuration 
would not be divergence free in the present setup. Fur- 
thermore, a field configuration facilitating wind ejection 
is expected to be generated automatically by the rotation 
of the disc ( Ouycd and Pudritz, 1997 ). Therefore, only the 
simple constant vertical magnetic field configuration has 
been implemented as initial condition. 



The magnetic flux density is specified by the parameter 
j3\ = P g as,i/-Pmag.i, at Tj = 1. In the numerical experiment, 
the density distribution is determined from Eq. H with 
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K = (7 - l)/ 7 - Using P gas = Kp-y and P raag - 5 2 /2 it 
follows that 



Bo = 



2K 



(11) 



£.3. Boundary conditions 



The disc boundary is located in the x = plane and is 
assumed to be a fixed base for the jet. As such, the ro- 
tational velocity, mass density and internal energy den- 
sity is kept constant at their initially prescribed values. 
Jet fluxes of mass, energy and momentum are assumed 
to have negligible impact on the disc surface layer. For 
instance, the disc is assumed to supply mass to the sur- 
face at the same rate as the jet carries mass away. This 
is implemented as a steady condition where the mass flux 
p x is symmetric at the lower boundary^] The energy flux 
at the disc surface has been implemented as either a zero 
x-derivative condition or a cold condition. The cold condi- 
tion is applied by assuming a jump in internal energy per 
unit mass (temperature) between disc and corona. The in- 
ternal energy density in the disc ghost zones are found by 
antisymmetrizing the a;- face centered values of T = e/p 
around (e /po)/£; 



T[b - i + 1/2] = 2(e /p )/e -T[b + i- 1/2], 



(12) 



with b = 4 and i = 1,2,3,4 so that the values in square 
brackets are the grid locations. The fixed internal energy 
density and mass density at the disc surface is denoted en 
and po respectively. The temperature jump is specified by 
the ratio £ = (T CO rona/Tdisc)o evaluated at the disc-corona 
interface. In effect, a cold inflow condition is specified by 
the internal energy density flux, h x = p x T. 

The inflow velocity is not specified, but allowed to 
evolve freely according to the forces acting on the disc 
surface. However, the a;-flux of ^-momentum (the xx- 
component of the Reynolds stress tensor) is not allowed to 
change the outflow velocity on the surface. Instead, ther- 
mal and magnetic forces are taken to be the dominant 
mechanisms for injecting gas into the corona. The pres- 
sure in the disc (ghost zones) and on the boundary is kept 
fixed at the values prescribing the initial pressure balance. 
Any pressure gradient at the disc surface is caused solely 
by deviations from the initial hydrostatic equilibrium oc- 
curring in the corona. 

For the evolution of the magnetic field the disc is taken 
to be perfectly conducting. In general, the magnetic field 
is only specified initially and is subsequently evolving ac- 
cording to the electric field. Keeping this in mind when 
imposing boundary conditions, violation of the divergence 
free condition of the magnetic field may be prevented in 

1 In the 1-D case, a symmetric or zero ^-derivative condition 
on the mass flux p x results in a steady state where dp/dt = 
—dpx/dx = at the boundary. 



a natural way. Particularly, one finds that the disc veloci- 
ties determine the electric field and thereby the evolution 
of the magnetic field at the lower boundary. The electric 
field on the disc surface is specified as 



E y — —u z Bq, 
E z = u y B , 



(13) 



where Bq is the constant vertical field specified initially by 
Eq. [TI]. The issues of boundary driving is discussed further 
in Sect. pTT . 



At the upper boundary, an extrapolation of the elec- 
tric field would not be stable and another approach has to 
be applied which is discussed in Sect. 2.3.2 . For the mass, 
momentum and energy density fluxes, a simple extrapo- 
lation into the ghost zones may be applied to allow the 
densities to evolve on the boundary. 

2.3.1. Boundary driving 

At the lower boundary the disc is kept rotating with a 
fixed angular velocity. This should in principle be sufficient 
to determine the evolution of the y- and ^-components of 
the magnetic field. It turns out, though, that the veloc- 
ity driving has to be specified in a somewhat non-local 
manner. This is because of the non-local nature of the 
numerical differential operators used in the code. The nu- 
merical scheme involves the three nearest points on each 
side along the direction of differentiation. Therefore, to 
avoid ripples when driving, the driving has to be passed 
smoothly to the active grid. To be specific, the velocity 
field just inside the active grid is determined by passing 
a third order polynomial through the boundary and the 
two neighboring zones further inside the grid, 



^[5] = gUi[4] 



6 Ui[6] - l;Ui[7], 



9 



9 



y,z, 



(14) 



where the lower boundary is located at zone index 4. The 
above chosen polynomial representation has been found 
in other but similar experiments to give the smoothest 



driving (Nordlund and Galsgaard, 1997a). It has not been 
subjected to tests in the present work. 

Care has to be taken also to avoid ripples caused by the 
periodicity of the y- and z-boundaries and by the shape 
of the box. To avoid shear at the box sides and reversed 
vorticity at the box corners, the disc velocity profile has to 
be terminated inside the periodic boundaries. As above, 
the non-local nature of the differential operators demands 
the decrease of disc rotation to span a few grid points. 
For the same reason, the velocity cutoff has to take place 
at some distance from the boundaries to ensure a region 
of practically zero rotation. To satisfy these conditions, a 
hyperbolic tangent function is used to cut off the velocity; 



fi(r) 



tanh 



1 



y,z. 



(15) 
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Fig. 1. Plots of the rotational velocity in the disc for different 
values of the softening parameter, l s . The solid, dotted, dashed 
and long dashed lines represent values of l s — (Keplerian), 
1/2, \/l/2 and 1 respectively. The plotted profile is given by, 



u(r) = uef where ug = r/(r + I 



2S3/4 



The cutoff function / 



terminates the rotation at the periodic boundary (r > 10) but 
does not affect the inner regions of the disc. 



Here Li is the size of the box in the i direction and s; is 
the distance that the midpoint of the hyperbolic cutoff is 
shifted from the boundary. The resulting velocity profile, 
Ui.cut = Uifi, is shown in Fig. |l] for L y = L z — 30 and 
s v = s z — 2.5 in a slice along y — and z > 0. 



2.3.2. The quasi-transmitting boundary 

A symmetric condition on the electric field does not allow 
tangential components of the magnetic field to evolve on 
the boundary. In particular, when the initial winding of 
the vertical magnetic field reaches the upper boundary as 
a toroidal Alfven wave, it would be reflected. To minimize 
reflection (reversal of the toroidal magnetic field) and pos- 
sible disruptive and/or decelerating effects on the flow, an 
upstream sensing of changes in the magnetic field is used 
to guide the evolution of the boundary field. This may 
appear somewhat artificial, but it corresponds roughly to 
apply values along an outgoing characteristic of a linear 
torsional Alfven wave. Such a wave propagates with the 
Alfven speed ma = B x /y/p. As a fair estimate, changes in 
the tangential field are taken to propagate with the (fixed) 
Alfven speed determined by 



(#o-*)' 



(16) 



evaluated at the upper boundary. 

The tangential magnetic field on the boundary is as- 
sumed to evolve towards some upstream value at a dis- 
tance Ax up , below the boundary. The travel time is es- 
timated from the Alfven speed; t up = Ax up /ua,x- The 



condition implemented on the electric field specifies the 
evolution of the tangential magnetic field as 

dBj = Afy 

8t t u r 



(17) 



The work done on the field may be controlled by the pa- 
rameter x, which specifies the convergence values for the 
magnetic field at the upper boundary; 



ABi.up = Bi[m x -4] -xBi 



i = y,z. (18) 



Here, n x specifies the index width of the sensing distance. 
Choosing X ~ 1 mimics in a simple way the work done 
on the magnetic field by the part of the jet outside the 
computational box. 

It must be emphasized, that the x parameter is to be 
regarded as a real physical parameter expressing to what 
extent the exterior is "braking" the field rotation. If, for 
example, the field is anchored to massive regions in the 
ambient medium outside the computational domain, the 
braking will cause the field to become inclined with respect 
to the outer boundary. Hence, the choice of x is n °t only 
a question of numerics, but reflects a real physical degree 
of freedom. 

2.4- Note on approach 

The model philosophy adopted here is different from pre- 
vious work, particularly in the applied boundary condi- 
tions. Instead of using large efforts on designing (artifi- 
cially) open outer boundary conditions, periodic bound- 
aries have been applied. The argument being, that it is 
inconceivable that real jets would not be perturbed by mo- 
tions in the surrounding medium anyway. In the present 
work, the outer boundaries as well as numerical noise pro- 
vide the perturbations triggering the possible instabilities 
we wish to investigate. Obviously, this approach cannot 
provide precise determinations of e.g. growth rates, but it 
is used to allow instabilities to develop in a natural way 
and provide indications on whether our understanding of 
jet collimation needs to be revised or not. 

3. Results 

The numerical experiments all show the same overall evo- 
lutionary sequence which may be summarized as follows: 

1. Winding and opening of the magnetic field lines. 

2. Relaxation of initial magnetic acceleration and spuri- 
ous magnetic reflections. 

3. Build up of collimated (but unsteady) outflow. 

Initially, the winding of the magnetic field propagates 
through the computational box as a torsional Alfven wave. 
The winding propagates with velocities expected to vary 
with distance approximately as Eq. [l^. Accordingly the 
winding propagates faster at large radii as noted also by 
Ouyed and Pudritz (|l997|) . 
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The upper boundary condition on the magnetic held 
allows (partially) for the Alfven pulse to be transmitted. 
Some reflections are unavoidable, but eventually the re- 
flections are transmitted through the boundaries and the 
system settles down. The relaxation of the minor Alfven 
pulse reflections concludes the two initial phases of the nu- 
merical experiments. At which evolutionary stage a (sig- 
nificant) collimated outflow is initiated in a real accreting 
system can only be speculated upon. In the present type of 
jet experiment, the initial phases do not reveal much rele- 
vant physics and the initial conditions as well as boundary 
conditions are probably not well suited for such investiga- 
tions. However, the setup eventually provides a collimated 
outflow with the expected qualities that we want to inves- 
tigate. 

A set of parameters have been introduced to truncate 
the physical problem and adopt it into a numerical setup. 
To probe the effect of these parameters, such as softening, 
box size and velocity cutoff, a range of experiments were 
carried out. The truncation of the problem and the effect 
of the related parameters have been presented above. In 
what follows, we will concentrate on a small number of 
experiments which have been conducted particularly to 
investigate stability issues. These experiments are listed 
in Table 0. 



run 


m x x m y x m z 


Lx X Ly X Lz 


ft 


t 


tin ax 


i 


98 x 71 x 71 


40 x 30 x 30 


5 


1 


500 


a 


98 x 71 x 71 


40 x 30 x 30 


10 


1 


500 


Hi 


98 x 71 x 71 


40 x 30 x 30 


5 


100 


500 


iv 


128 x 85 x 85 


60 x 40 x 40 


10 


1 


400 


V 


200 x 101 x 101 


200 x 25 x 25 


10 


1 


260 



Table 2. List of experiments for the stability study. All ex- 
periments use the softening parameter, l s — v/l/2, and a poly- 
tropic index, 7 = 5/3. 



3.1. Jet dynamics 

By examining the rates of change of energies, the action of 
various forces and issues of stationarity may be analyzed. 
The rate of change of total energy in a volume of gas is 
controlled by the rate of energy transport in and out of 



the volume, i.e. the net flux, AF = F(x\ owcr ) — F(a 



r), 



the rate of energy conversion through work, W, and dissi- 
pation, Q. Consequently, the rate of change of magnetic, 
A4, kinetic, K,, thermal, T, and gravitational energy, Q, in 
a volume may be expressed as follows, 



dM 



AF maR + W Lo 



— AFkin + AFvisc + Wgrav- 



gas 



(19) 



' WLorentz — Qvisc 7 (20) 



dT 

-rrr — AF ent h 
at 

^-AF - 
dt ~ grav 



5 Joule 1 



grav ■ 



(21) 



(22) 



The rate of total work done by gravity in the box, W grav , 
the rate of work done against gas pressure gradients, W gas , 
and the rate of work done against the Lorentz force, 
WLorentz, are given by 



grav 



u ■ pV<S> dV, 



Wg as = / u-WW, 



Wlo 



- / u • (J x B) dV. 

Jv 



(23) 



(24) 



(25) 



Here dV = dx dy dz is an element of the volume, V = 
flower, £ U pp Cr ] x [-L y /2,L y /2] x [-L z /2,L z /2] . The en- 
ergy fluxes through a cross section, x, are given by 



F mag (x) = / (E y B z - E z B y ) dS x , 
Js 



firing) 



1 



pu dS x , 



F eat h(x) = / pu x jedS x 
Js 



-^grav (^) 



- grav 

where dS x 
face normal 



pu x $dS x , 



(26) 
(27) 
(28) 
(29) 



e dy dz represents the element of a sur- 
to the ^-direction and the surface, S = 
[— Ly/2, L y /2] x [— L z /2, L z /2], is an entire cross section 
of the box. 

To demonstrate the quasi-stationary state of the flow 
and the general energy conversion mechanisms, the time 
averaged energy fluxes (Eqs. p6]-|29]) are plotted in Fig. || 
as functions of height above the disc. The fluxes are ob- 
tained as differences between the curves as indicated by 
the arrows in the plots. 

For both experiments it is seen that more magnetic 
energy is transported into the volume than is carried out. 
The magnetic energy does not increase accordingly, since 
the difference between magnetic energy flux in and out 
of the volume is balanced on average by Lorentz work 
and Joule dissipation. The conversion of magnetic energy 
flux into kinetic energy flux is dominant close to the disc, 
whereas magnetic energy flux is primarily converted into 
thermal energy flux in the upper two thirds of the box for 
both the experiments. Just above the disc surface, a (hy- 
drostatic) balance is seen to be maintained between grav- 
itational energy flux and the sum of enthalpy and kinetic 
energy flux. This corresponds to a situation where gravity 
is balanced by gas pressure and x-advection of momentum. 
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.^grav 
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Fig. 2. Time averaged total energy fluxes as functions of height for run ^] (left panel) and run |^ (right panel). The dotted line 
is magnetic flux, dashed line is the sum of kinetic and magnetic fluxes, the dot-dashed line is the sum of enthalpy, kinetic and 
magnetic fluxes. The solid line is the total energy flux, which is obtained by adding the (negative) flux of gravitational energy 
to the sum of the other fluxes (dot-dashed line). 



The two experiments differ significantly in the amount of 
magnetic energy flux converted into thermal energy flux 
as a consequence of the difference in the magnetic energy 
reserve at hand. 

Details of the energy conversion processes as a func- 
tion of time are shown in Fig. |^. After a sharp increase 
of kinetic energy, caused by the initial Alfven pulse, the 
rate of change of kinetic energy oscillates as does the rate 
of work done by the Lorentz force. Both the efflux of ki- 
netic energy and the work done by gravity are fairly con- 
stant, but large variations in the work done by gas pres- 
sure are seen for run pi The characteristics of the events 
at the "pressure work peaks" (t w 270, 400, 450) for run § 
are all indications of major, dynamically important mag- 
netic reconnection events. Such events relax the wound up 
magnetic field configuration, whereby Lorentz work is re- 
duced as there is no longer the same amount of azimuthal 
magnetic field for driving (and pinching) the flow. Conse- 
quently, the rate of change of kinetic energy drops. Before 
this happens, the thermal energy has been building up for 
some time (positive rate of change of thermal energy) and 
much of this energy is now released by pressure work. This 
promptly makes the rate of change of kinetic energy pos- 
itive. As the magnetic pinching is significantly reduced, 
one possible effect of the pressure work "explosions" is 
to generate filaments and disruptions of the jet into the 
ambient medium. The magnetic pressure of the surround- 
ing vertical magnetic field eventually halts this (irregular) 
radial expansion (cf. Fig. ^). The events at t s» 400 and 
t w 450 are followed by an increase in the rate of Lorentz 
work which supports this picture. As a consequence of the 
work done by gas pressure, and the increased enthalpy 
flux out of the volume following the reconnection events, 
the rate of change of thermal energy becomes negative. 
The slight decrease in net Poynting flux is in line with 



a constant Poynting flux at the lower boundary but a 
decreased Poynting flux out of the volume. The sudden 
deficit in Poynting flux through the upper boundary is a 
consequence of azimuthal field relaxation in the box, caus- 
ing the tangential field components to decrease at exit. 

Comparing the two experiments run |^ and run ^| the 
weak field experiment (run appears much less violent 
with respect to magnetic reconnection events. The most 
prominent events to be identified occur relatively late, at 
t ~ 320 and t w 400. Apparently these events do not cause 
significant changes in the jet dynamics, as no abrupt peaks 
in the pressure work are detected. Instead, another promi- 
nent feature of the jet dynamics may be identified, namely 
the oscillatory pattern in the rate of Lorentz work and ki- 
netic energy flux. Such harmonic oscillations (in the radial 
flow) between an inner ram pressure barrier (centrifugal 
barrier) and an outer magnetic pressure barrier have been 
predicted analytically ( |5auty and Tsinganos, 1994 ) and 
shown in axisymmetric experiments to steepen into fast 



magneto-sonic shocks (Ouyed and Pudritz, 1997). The os- 
cillations are most clearly present until approximately the 
time where the first prominent reconnection event takes 
place and complicates the flow pattern. 

3.1.1. Forces and flow features 

A look at the forces reveals that close to the central object 
the gas is launched from the disc surface by a thermal pres- 
sure gradient between disc and corona. Just above the disc 
the acceleration process may be regarded as either mag- 
netic and centrifugal. Close to the rotation axis the mag- 
netic point of view may conveniently be adopted as the 
magnetic field is here highly wound up and a vertical gra- 
dient of the field is the main acceleration mechanism. At 
larger radial distances, the magnetic field is not too wound 




Fig. 3. The rate of change of magnetic, kinetic and thermal energy in the box between x = 1.1 and x = 38.8 are graphed 
for run | (left panels) and run ^ (right panels). The top panels show rate of change of magnetic energy (dotted line), the net 
Poynting flux into the volume (solid line) and the rate of work done by the Lorentz force (dot-dashed line). The middle panels 
show rate of change of kinetic energy (dotted line), net kinetic flux into the volume (solid line) and the rate of work against 
gravity done by carrying matter through the volume (dot-dashed line) . The lower panels show rate of change of thermal energy 
in the volume (dotted line), net enthalpy flux (solid line) and the rate of work done by gas pressure (dot-dashed line). The 
dashed line marks zero rate of change of energy. Please note that until t — 150 only 10 snap shots were saved (at time intervals 
At = 15). 
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up for the gas to be flung more radially outwards. Here the 
acceleration mechanism bares the same characteristics as 
in the centrifugally driven wind scenario. The central jet 
region appears very well collimated from the start, as the 
magnetic acceleration process for this region is inherently 
vertical and causes little radial flow. The fraction of the jet 
flung more radially outwards is eventually collimated by 
the magnetic force. Fig. || shows the collimating magnetic 
force decomposed into toroidal and radial magnetic pres- 
sure. The ambient vertical magnetic field has generally a 
small collimating effect, but when magnetic reconnection 
events relaxes the wound up field structure the vertical 
field takes over and prevents locally and for a period of 
time the flow from expanding into the ambient medium. 

Radial magnetic force 
0.06 1 1 1 1 1 l 1 1 1 1 i 1 1 1 1 1 1 1 i 1 1 1 1 i 1 1 1 1 



0.04 - 



0.02 - 



-0.02 - 



-0.04 - 




-0.06 I i i i i I i i i i I i i i i I i i i i I i i i i I i i i i 

-10 -5 5 10 

y 

Fig. 4. The collimating magnetic forces at x ~ 4 for run ^ at 
t = 200. The dashed line is the radial (here the y-direction) 
pressure gradient of the toroidal magnetic field (here, B z ) and 
the dotted line is the radial pressure gradient of the vertical 
magnetic field. 







N 



-5 



y 

Fig. 5. Complexity of the flow illustrated in a cross section at 
x « 20 for |i| at t = 270. The dashed line contours mark zero 
vertical velocity, dash-dot-dot-dot line contours mark the sonic 
points and the solid line contour marks the Alfven surface. 

To illustrate the complexity of the resulting flow mo- 
tions, Fig. H shows velocity vectors in a cross section of the 
flow. Contours indicating the Alfven surface, super-sonic 
regions and regions of back-flow are drawn. The super- 
sonic and super- Alfvenic jet beam is surrounded by a su- 
personic shell at larger radii. Furthermore, a prominent 
back-flow region is noted just outside the jet beam in the 
upper right quadrant of the plot. 



Velocity pattern at x=20.2696 




The mean vertical velocity increases approximately lin- 
early with height and the terminal velocities obtained in 
the jet at the upper boundary are of the order of the Ke- 
pler velocity at the inner disc radius. However, there is 
some dependence of the mean terminal velocities on the 
strength of the initial field as expected from the predic- 
tions by Michel (Michel, 1969; see also Kudoh et al., 1997). 
The time averaged terminal velocities are consistent with a 

2/3 

power law dependency, Uoo cx B ' , in the relatively early 
quiescent stages of the jet evolution. At later stages, dis- 
ruptive events become dynamically important and the as- 
sumptions for the predicted terminal velocity dependency, 
e.g. that the gas dynamics is governed purely by magnetic 
effects, do not hold. The more violent disruptive events 
in the strong field experiment redirect a relatively larger 
fraction of the vertical flow into transverse gas motions. 
Eventually, the mean vertical velocity of the strong field 
experiment becomes less than the mean vertical velocity 
of the weak field experiment. 



3.2. Jet stability 

Though potentially disruptive events may be identified in 
the dynamics, these events are found only to generate fil- 
aments and cause non-destructive distortions of the jet 
beam. To monitor the onset and evolution of instabilities 
further, the magnetic field topology is investigated. Phc- 
nomenological investigations of the evolution of the field 
topology are carried out by visualizing isosurfaces of mag- 
netic flux density. Fig. |?] shows the evolution of an iso- 
surface of the magnetic flux density. The magnetic field is 
seen to develop a complex and highly non-axisymmetric 
structure, which relaxes periodically allowing the field to 
return to a less complex, more nearly cylindrical topology. 
In each such build- up/relaxation cycle several modes may 
be identified qualitatively such as the sausage instability 
and the elliptical (\m\ = 2) modes. 

The total viscous and magnetic dissipation as a func- 
tion of time are shown in Fig. || for a test volume. In 
the initial phase the wound up field structure is build- 
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Fig. 7. Isosurfaces of magnetic flux density for run |. The snap shots are taken at t = 270, 274, 278, 282, 286, 290 starting at 
the upper left panel and ending at the lower right. The middle section of the flux tube is located off center and partly outside 
the test volume at first, but the non-axisymmetric distortion relaxes and the tube swings (clockwise) into a more axisymmetric 
configuration with a tightly wound spiral structure close to the disc surface. For clarity, only a narrow region of the computational 
domain centered around the a;-axis is shown. 



ing up and the magnetic dissipation is seen to increase Fig. [7] is seen to be connected to a significant increase in 

sharply until t ps 100. From the magnetic dissipation, sev- the magnetic dissipation. A slight increase is seen also in 

eral minor and a couple of major magnetic reconnection the viscous dissipation, indicating that a small fraction of 

events are identified to occur in this region for run || (up- the energy released by field re-organization is transferred 

per panel). The evolution of the magnetic field shown in into kinetic energy and subsequently dissipated. 
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/lagnetic field orientation at x=20.2696 




Fig. 6. The magnetic field orientation for run a at t — 270 in 
a cross section at i « 20. The contour line marks the Alfven 
surface. 



Dissipation in test volume for ^=5 




100 200 300 400 



Dissipation in test volume for /? ; =1( 




Fig. 8. Total dissipation in a test volume for run |/| (upper 
panel) and run ^ (lower panel) . The dissipation is integrated 
in a central volume of the box between x = [0.7,38.1] and 
y,z = [—4.7,4.7] and is normalized by the maximum value of 
the magnetic dissipation in each experiment. The solid lines 
are magnetic dissipation and the dashed lines are viscous dis- 
sipation. 



tant consequences for the jet stability. The significance of 
the magnetic diffusion varies in correspondence to e.g. the 
characteristic flow velocity as. Specifically, one expects in 
the Pi = 10 case (run ^ the characteristic (vertical) jet 
velocity to be less than in the $ = 5 case (run |) as noted 
in Sect. |3.1.1| . Accordingly, the magnetic Reynolds num- 
ber, ReM — UL/rj, will be smaller in the weak field exper- 
iment and magnetic diffusion relatively more important. 
The magnetic diffusion is better capable of counter bal- 
ancing the continual field winding and the build-up phase 
of the azimuthal field is prolonged. For the same reason 
the relaxation events themselves appear less violent, as 
seen by comparing the panels of Fig. |^, and will cause less 
destructive, more localized distortions of the flow. 

The damping of the spiraling jet motion at heights 
where the jet is in general super- Alfvenic is related to 
what appears as field "unwinding" just outside of the 
super-Alfvenic central beam. Fig. |] shows the orientation 
of the magnetic field in a cross section at x ss 20 and il- 
lustrates the unwinding. The jet beam swings clockwise in 
the plot, and the change in field orientation is seen partic- 
ularly on the front of the Alfven surface in the direction 
of the (clockwise) spiral motion. 

The unwinding of the field comes about as the spi- 
ral motion of the jet beam forces the wound up field into 
the ambient almost vertical field. The wound up field is 
oriented in a direction practically perpendicular to the 
ambient vertical field and reconnection occurs when these 
oppositely directed field lines are forced into a small re- 
gion. This is seen as sheets of Joule dissipation marking 
the regions of magnetic diffusion at the jet flanks in the 
left panel of Fig. || The right panel of Fig. ^ shows such a 
sheet in more detail and the orientation of the field lines 
reconnecting. These magnetic reconnections, occurring in 
a "cocoon" surrounding the central jet beam, are the rea- 
son for the observed back-flow in Fig. || and the field un- 
winding in Fig. [| More specifically, the upper right region 
of back-flow in Fig. ^| corresponds to where the recon- 
nection events are expected to catapult gas backwards in 
this scenario. The gas at the "nose" of the Alfven surface 
(at y Rj 2, z rs —1) is accelerated oppositely, i.e. forward 
in the rotation direction and upwards. Also seen in left 
panel of Fig. |^ is a central region of significant magnetic 
dissipation. This is the region where reconnection occurs 
as a consequence of field interlocking resembling the sce- 
nario observed in experiments c oncerning coronal heating 
( Galsgaard and Nordlund, 1996 ). The right panel of Fig. [| 
confirms this interpretation. 



The dissipation of the magnetic field does not match 
the Poynting flux into the test volume in general. The 
surplus of magnetic energy is partly used for driving the 
jet flow, as demonstrated in Sect. 3.1, and partly to build 
up the winding between relaxation events. In addition, 
there is a continuum of "background" magnetic dissipa- 
tion between the relaxation events which may have impor- 



4. Discussion 

Obviously, parameter space has not been probed in all 
detail in the work presented here. Issues which clearly 
need more attention in future work include the magnetic 
configuration, thermal conditions and mass loading. Prob- 
ing other magnetic configurations would involve a signifi- 
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Fig. 9. The left hand side panel shows isosurface of magnetic dissipation. The dissipation is seen to occur predominantly in 
two regions; in the central highly wound field region close to the rotation axis and in a "cocoon" at the jet flanks. The right 
hand side panel shows zoom-in on magnetic dissipation isosurfaces and inner wound up magnetic field field lines pressed towards 
the outer more vertical field. Both snap shots are from run m at t — 300. 



cant change of the numerical code, whereas the inclusion 
of detailed cooling expressions would be straightforward 
( Rognvaldsson, 1999 ; Rognvaldsson et al. , 2002 ) . Magnetic 
field configurations which fan away from the jet axis in a 
more dipole like fashion may provide less stabilization high 
above the disc and far from the rotation axis. However, 
close to the central parts of the disc the field lines will, in 
a pote ntial configuration as suggested by Cao and Spruit 
(1994), bend towards the jet axis as a function of height 
and not appear much different from the vertical configu- 
ration used in this work. The stabilization and collimation 
provided by the poloidal field may be reduced if the mag- 
netic field is produced in a dynamo active disc. In such a 



scenario, Brandenburg et al. (2001) have found that signif- 
icant amounts of toroidal magnetic field are transported 
into the corona from the outer parts of the disc. This could 
seriously influence the stability. 

Theoretically, the stability of twisted magnetic flux 
tubes is expected to depend on the diameter to length 
ratio (aspect ratio), such that relatively long tubes are 
most unstable. This is strictly valid only in the ideal mag- 
netohydrodynamics case. In practice, the rate of magnetic 
dissipation tends to increase for more narrow tubes, which 



reduces the importance of the aspect ratio (Galsgaard and 



Nordlund, 1997). In the work presented here, no system- 



atic change was found in the appearance of the jet when 
the aspect ratio of the experiment was changed (e.g. com- 
paring runs |^, ^ and 0) . 



As a consequence of magnet ic d issipation occurring 
primarily in two regions (cf. Sect. 3^2) the jet flow consists 
of a hot central beam, with temperature of the order of 
the virial temperature at the inner disc radius, and a hot 
cocoon. Hot jet flanks has been observed in YSO's and 
are proposed to be generated by shocks associated with 
time variability and entrainment in the flow (e.g. HH47, 
Hartigan et al., 1993). The results presented in Sect. ^2 
suggest that magnetic dissipation may be a major addi- 
tional mechanism for heating in this region. 

Preliminary results (run indicate that relatively 
cold disc outflow results in less distorted, more well- 
defined flows and relatively more dense jets. A similar re- 
sult has been reported by Hardee et al. (1997) who found 
dense jets to maintain a high-speed spine which prevents 
disruption of the internal jet structure though large helical 
and elliptical distortions were present. The experiment, as 
it is, is already suitable for addressing further which con- 
sequences the temperature of the disc outflow may have 
on the jet dynamics. However, further investigations with 
special emphasis on the transsonic region close to the disc 
surface are desirable to investigate the issues of mass load- 
ing and disc-jet interactions in general. 



5. Summary and conclusions 

A high order numerical scheme has successfully been 
adopted and a suitable mesh refinement specified. Ini- 
tial conditions resembling previous axisymmetric numeri- 
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cal experiments have been chosen to ease comparison. The 
polytropic equation of state is only used initially, to pre- 
scribe the initial density distribution of the corona. The 
most important features that differ from previous jet ex- 
periments are: 

— The model is three dimensional rather than axisym- 
metric. Due to the periodic boundaries, a cutoff of the 
disc at large radii is applied. 

— The thermal energy equation is solved, with self- 
consistently computed heating by viscous and Joule 
dissipation. 

— "Free" mass outflow from the disc, i.e. the mass flux 
is allowed to adjust self-consistently to the forces near 
the disc surface. 

— Parameterized Poynting flux through the upper 
boundary, representing external work. 

— The experiments have been evolved far beyond the ini- 
tial transient, and display a quasistationary behavior, 
as evidenced for example by the nearly constant total 
energy flux. 

The jet dynamics has been investigated by analyzing 
the mechanisms of energy conversion. The rotational en- 
ergy of the disc is carried by the magnetic field into the 
corona and is first predominantly converted into kinetic 
energy. In the upper two thirds of the computational do- 
main the magnetic energy is predominantly converted into 
thermal energy. General features predicted by steady state 
theory and axisymmetric numerical experiments, such as 
knot generation and terminal velocity dependency on the 
magnetic field strength have been confirmed in the rela- 
tively early and quiescent stages of the experiments. At 
later stages the flow becomes quite unsteady as instabil- 
ities set in, but no serious disruption of the flow occurs. 
The jet stability is found to be influenced by the magnetic 
dissipation — this has not previously been investigated 
in the context of jet flows. The main findings concerning 
magnetic dissipation are the following: 

— The heating by magnetic dissipation is significant, and 
leads to jet temperatures of the order of the virial tem- 
perature of the innermost Kepler orbit. 

— Magnetic reconnection occurs primarily in two regions: 
in a central region of the jet due to field interlocking 
and in a jet cocoon due to the spiraling motion of the 
jet beam forcing the wound up field into the ambient 
vertical field. 

— Magnetic dissipation helps to prevent critical kinking, 
and the jet is able to sustain a quasistationary flow, 
with only limited excursions. 

— Reconnection events are seen to result in mass ejection 
into the ambient medium and cause filament structures 
in the jet beam. 

— Reconnection events are seen to significantly influence 
the dynamics at the jet flanks where deceleration and 
even back-flow can be found quite close to the central 
super- Alfvenic beam. 



— Heating in the region just above the disc is likely to 
have a significant effect on the mass loading. 
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